//
// Distributions.cxx
//
// Created by TB on 11/30/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <iostream>

#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>

#include "tppmc/CollinsFragmentation.h"
#include "tppmc/Configuration.h"
#include "tppmc/ConstantFragmentation.h"
#include "tppmc/Distributions.h"
#include "tppmc/DssFragmentation.h"
#include "tppmc/Dssv.h"
#include "tppmc/CollinsKt.h"
#include "tppmc/Grv98Pdf.h"
#include "tppmc/sivers/Torino2009.h"
#include "tppmc/sivers/Positivity.h"
#include "tppmc/SofferBound.h"
#include "tppmc/SofferBoundProkudin.h"
#include "tppmc/TransversityProkudin.h"

#ifdef USING_LHAPDF
#include "tppmc/LhapdfDistribution.h"
#endif

namespace tppmc {

typedef boost::shared_ptr<PartonDistribution> PartonPtr;
typedef boost::shared_ptr<FragmentationFunction> FragPtr;

Distributions::Distributions() {
   // Assign PDF set. No option currently.
   mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
   // Assign Sivers set. Options are Anselmino extraction or
   // positivity bound (needs to be implemented).
   mSivers = PartonPtr(new sivers::Positivity());
   // Assign transversity set. Options are Prokudin extraction
   // or Soffer bound.
   mTransversity =
      PartonPtr(new SofferBoundProkudin(SofferBoundProkudin::LO));
   // Assign fragmentation function. Options are constant fraction
   // of parton momentum, DSS, DSS + Collins, KKP and AKK (KKP
   // and AKK need to be implemented).
   mFragmentation = FragPtr(new ConstantFragmentation(1.));
}

Distributions::Distributions(const Configuration& config) {
#ifdef USING_LHAPDF
   if(config.LhapdfSet() > 0) {
      mUnpolarised = PartonPtr(new LhapdfDistribution(config.LhapdfSet(),
                                                      config.LhapdfMember()));
   } // if
   else {
      mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
   } // else
#else
   mUnpolarised = PartonPtr(new Grv98Pdf(Grv98Pdf::LO));
#endif
   // Assign Sivers set. Only choice is Anselmino extraction or-
   // positivity bound based on the unpolarised distribution.
   switch(config.SiversType()) {
      case 1: // Anselmino et al. 2009, normal (not positivity bound)
         mSivers = PartonPtr(new sivers::Torino2009(false));
         break;
      case 2: // Positivity bound based on unpolarised PDF
         mSivers = PartonPtr(new sivers::Positivity(mUnpolarised));
         break;
      default: // No sivers
         mSivers = PartonPtr(new NullDistribution);
         break;
   } // switch
   
   PartonPtr dssv(new Dssv);
   switch(config.TransversityType()) {
      case 1: // Anselmino et al., Phys. Rev. D 75, 054032 (2007)
         // Select either leading or next-to-leading order
         if(config.Order() == 0) {
            mTransversity = PartonPtr(
               new TransversityProkudin(SofferBoundProkudin::LO));
         } // if
         else {
            mTransversity = PartonPtr(
               new TransversityProkudin(SofferBoundProkudin::NLO));
         } // else
         break;
      case 2: // Soffer bound from GRV98 and GRSV
         // Select either leading or next-to-leading order
         if(config.Order() == 0) {
            mTransversity = PartonPtr(
               new SofferBoundProkudin(SofferBoundProkudin::LO));
         } // if
         else {
            mTransversity = PartonPtr(
               new SofferBoundProkudin(SofferBoundProkudin::NLO));
         } // else
         break;
      case 3: // Soffer bound from chosen unpolarised PDF and DSSV
         mTransversity = PartonPtr(new transversity::SofferBound(
                                   mUnpolarised, dssv, false));
         break;
      default: // No transversity
         mTransversity = PartonPtr(new NullDistribution);
         break;
   } // switch
   
   // Retained support for --quark_transv_soffer for backward compatibility.
   // Assign transversity set. Options are Prokudin extraction
   // or Soffer bound.
   if(config.UseTransversitySofferBound()) {
      mTransversity = PartonPtr(new transversity::SofferBound(
                                mUnpolarised, dssv, false));
   } // if

   // Assign fragmentation function.
   // If fragmentation is not requested hadrons are generated with a
   // constant fraction of the parton momentum.
   // Options for fragmentation are DSS, DSS + Collins, KKP and AKK (KKP
   // and AKK need to be implemented).
   switch(config.FragmentationFunction()) {
      case 0: // No fragmentation
         mFragmentation = FragPtr(new NullFragmentation());
         break;
      case 1: // DSS
         mFragmentation = FragPtr(new DssFragmentation(config.Order(),
                                                       config.MeanJt2()));
         break;
      case 2: // Constant fragmentation
         mFragmentation = FragPtr(new ConstantFragmentation(1.));
         break;
      case 4: // DSS + Collins
         mFragmentation = FragPtr(new CollinsFragmentation);
         (dynamic_cast<CollinsFragmentation*>(mFragmentation.get()))
         ->SetPositivityBound(config.UseCollinsPositivityBound());
         break;
      case 5: // DSS+Collins with kT-dependent positivity bound
         mFragmentation = FragPtr(new collins::PositivityKt(FragPtr(
            new DssFragmentation(config.Order(), config.MeanJt2())),
            config.MeanJt2(), config.Hadrons()));
         break;
      default:
         mFragmentation = FragPtr(new ConstantFragmentation(1.));
         break;
   } // switch
   mUnpolarised->Initialise();
   mSivers->Initialise();
   mTransversity->Initialise();
   mFragmentation->Initialise();
}

Distributions::~Distributions() {
}

PartonPtr Distributions::Unpolarised() const {
   return mUnpolarised;
}

PartonPtr Distributions::Sivers() const {
   return mSivers;
}

PartonPtr Distributions::Transversity() const {
   return mTransversity;
}

FragPtr Distributions::Fragmentation() const {
   return mFragmentation;
}

std::list<PartonPtr> Distributions::Pdfs() const {
   std::list<PartonPtr> pdfs;
   if(mUnpolarised) {
      pdfs.push_back(mUnpolarised);
   } // if
   if(mSivers) {
      pdfs.push_back(mSivers);
   } // if
   if(mTransversity) {
      pdfs.push_back(mTransversity);
   } // if
   return pdfs;
}

void Distributions::Speak() const {
   std::cout << "Unpolarised:   " << Unpolarised()->Name()   << std::endl;
   std::cout << "Sivers:        " << Sivers()->Name()        << std::endl;
   std::cout << "Transversity:  " << Transversity()->Name()  << std::endl;
   std::cout << "Fragmentation: " << Fragmentation()->Name() << std::endl;
}

} // namespace tppmc
